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Abstract. We develop a recently proposed importance-sampling Monte Carlo algorithm for 
sampling rare events and quenched variables in random disordered systems. We apply it to a 
two dimensional bond-diluted Ising model and study the Griffiths singularity which is considered 
to be due to the existence of rare large clusters. It is found that the distribution of the inverse 
susceptibility has an exponential tail down to the origin which is considered the consequence of 
the Griffiths singularity. 



1. Introduction 

Monte Carlo (MC) methods are general computational techniques for sampling from a given 
probability distribution and for estimating expectation values under the distribution. They have 
been used in a wide range of physics pQ and statistical science [2]. Most of the MC methods are 
based on the Metropolis importance-sampling strategy, in which a Markov chain is constructed 
so that its stable distribution converges to a required distribution. Some improvements on 
the Metropolis MC algorithm have been made in order to accelerate the sampling process. 
One of them is extended ensemble MC methods [3] which includes multicanonical method [1], 
simulated tempering [5] and exchange MC method [6] (parallel tempering). These methods have 
been considered to be quite useful for studying complex systems such as spin glasses [71 151 l9l [10] . 
Among them, the multicanonical method successfully applied to statistical-mechanical models 
where rare events play an important role in the physical nature. Such examples of rare events 
are a mixed phase at a first-order phase transition temperature of the two-dimensional 10-state 
Potts model [I] and tails of an order-parameter distribution of the two-dimensional Ising model 

EH- 

Another class of rare events, which has not extensively studied yet, is found in quenched 
disordered systems when a distribution of a physical quantity over the quenched disorder is taken 
into account. A typical example is the rare-event tails of the ground-state-energy distribution of 
quenched disordered system that has been the subject of recent theoretical studies [HI [13]. It is, 
however, difficult for a conventional simple-sampling method to evaluate the tails of distribution 
precisely. Recently Korner, Katzgraber and Hartmann [16] have proposed an importance- 
sampling MC algorithm in order to evaluate the tails of the ground-state-energy distribution 



with high precision. Using the algorithm they could evaluate the tails of the distribution of a 
spin-glass model up to 18 orders of magnitude. 

A Griffiths singularity [Hj of random spin models is also caused by the rare events, namely 
the existence of arbitrary large clusters. Although a dynamical aspect of the Griffiths singularity 
has been confirmed by numerical simulations, there is, to our knowledge, no numerical evidence 
of the singularity in static observables and the function form that the Griffiths singularity takes 
has not been well established. It is the main purpose of the present work to detect numerically 
the Griffiths singularity. We develop the importance-sampling algorithm mentioned above and 
apply it to a finite-temperature calculation of the susceptibility distribution of a bond-diluted 
Ising model which is expected to exhibit the Griffiths singularity. 

The outline of this paper is as follows. In section [2] we explain the Monte Carlo algorithm 
for efficiently sampling the quenched variables of the disordered systems. In section [3] we briefly 
review the Griffiths singularity in randomly bond-diluted Ising model. In section 0] we study the 
Griffiths singularity by applying the method described in section [2j In section [5] we finally draw 
our conclusions and perspectives. 

2. Monte Carlo algorithms for quenched variables 

A quenched disordered system is defined by a Hamiltonian iT(S|J), where S denotes a 
configuration of the system and J a set of quenched variables generated by the probability 
distribution -P(J). In spin-glass models, for instance, S and J correspond to the spin variables 
and the interaction bonds, respectively. One needs to take double averages for evaluating a 
physical quantity v4(S,J) in the quenched disordered system. Firstly, a thermal average (A) (J) 
at temperature T = 1/(5 is calculated by tracing over the variable S for a fixed set of J, which 
is expressed as 

(A)(J) = ^jyTr s exp(-^(S|J)), (1) 

where the normalization factor Z(3) is the partition function. The other average is to take an 
average (A) (J) over the distribution P(3), 

[(A)] = J d3P(3)(A)(3). (2) 

In general, the distribution of the quantity A is defined by 

P(A) = J d3P(3)5 (A - (A)(3)) , (3) 

where 5 is the Dirac delta function. The aim of this work is to estimate accurately the 
distribution P(A), particularly the tails of P(A), for any physical quantity A we require at 
finite temperatures. 

2.1. Simple- sampling algorithm 

A standard method for estimating the distribution P(A) in equation ([3]) is to use a simple 
sampling MC algorithm, where P(A) is estimated as 

-, M 

P(A) = -Y, S ( A -( A K j{i) )) ( 4 ) 

with M being the number of the set of J generated independently from the given distribution 
P(3). This is only histogram accumulating the thermal average (^4) for each J. It is noted that 
in most quenched disordered systems, the thermal average is a difficult task at low temperatures 



because it takes a huge amount of time to equilibrate. The so-called extended ensemble MC 
methods [3] such as the multicanonical method [I] and the exchange MC method [6 J (parallel 
tempering) have partially overcome the difficulty of slow relaxation and/or have significantly 
reduced equilibration time. This, together with recent progress of computer power, allows us 
to take a large number of average over the quenched variables. A typical example of the total 
number of samples M would be of order of 10 3 in recent MC simulations of Ising spin glasses 
[TU]. Such a number is sufficient for evaluating the typical values, the averaged physical 
quantities, in most simulations. It is, however, still difficult for the simple-sampling algorithm 
to estimate rare cases, i.e., the tails of the distribution of equation ([3]). Particularly, it is 
impossible in principle to obtain the distribution whose value is less than 1/M. 

2.2. Importance- sampling algorithm for quenched variables 

In this work, we discuss an efficient MC method for sampling the quenched variables. A 
pioneering work was given by Hartmann [15], who studied a large deviation function of a 
sequence alignment. In the literature of physics, an importance-sampling MC algorithm for 
the quenched variables was proposed by Korner, Katzgraber and Hartmann and applied to the 
measurement of the ground-state-energy distribution in the Sherrington-Kirkpatrick (SK) model 
[T6] . Subsequently, the method was used for estimating the ground-state-energy distribution for 
the directed polymer in a random medium [17J. An application have been proposed to the 
problem of estimating the bit-error distribution of error-correcting codes [18] . 

An importance-sampling has been proposed for efficiently sampling the quenched variables J 
|15^I16]. The main idea for the importance-sampling is to enhance a probability for finding a set 
of J which gives a rare event for (.A) (J) in the distribution P(A) in Markov-chain MC simulations. 
In reference [16], a guiding function P(A) is introduced as a guess for the true distribution and 
the quenched variables are sampled from the inverse of the guiding function, 1/P{A) using the 
importance-sampling technique. When the guiding function P(A) approximates well the true 
distribution P(A) which is unknown a priori, the probability for visiting a set of J is nearly 
independent of the thermal average (A)(3) for a given J. This means that a histogram of the 
quantity A is nearly flat. The idea behind this algorithm is the same as the multicanonical 
method [1] where a guiding function as a function of the energy is introduced, instead of the 
standard Boltzmann weight, so that a resulting energy histogram becomes flat. 

The remaining problem for performing MC simulations is how to guess the guiding function 
P(A). In [16\ I17j . the guiding function is assumed to be a modified Gumbel distribution, in 
which a set of parameters is determined by fitting the data obtained by a preliminary simple- 
sampling run. The method does work quite well in the applications to the ground-state-energy 
distribution of the SK spin glass [16] and the directed polymer |17| . However it might rely 
on if the true distribution is expressed approximately by a specific chosen analytic function 
such as a (modified) Gumbel function. In this work, we use rather a learning algorithm of the 
multicanonical method, which is a recursion scheme proposed by Wang and Landau [191 120] . to 
obtain the guiding function with no assumption of the function form. 

The importance-sampling MC algorithm used in this work is defined by the following steps. 
Some initial estimate is made for the guiding function P(A); i.e., P(A) is constant or the 
resultant histogram by a simple-sampling run. An initial configuration of J is randomly chosen 
for a set of J according to the given distribution P(3). A modification factor / for the guiding 
function P(A) is initially set to a sufficient large value, e.g., f = e 1 ~ 2.718281 • • • [El [20] . 

(i) From the current i-th configuration jW, a new candidate J' is produced by replacing a 
subset of chosen at random with new values generated by the given distribution P(3). 

(ii) For the new candidate J', the thermal average (A) (J') is calculated by some numerical 
method that depends on the system we consider. We call this procedure inner loop in the 



sense that one has to take an average over the variable S for a given J, while the update 
scheme for J discussed here is referred to as outer loop. We shall discuss implementation 
of the inner loop later. 

(iii) The new candidate J' is accepted with probability 

P minf P((A)(jW) ll (5) 

using the current guiding function P(A) and the (i + l)-th configuration set to J'. 

When the new candidate is rejected the current one has to be counted again. 

(iv) The guiding function P(A) with A = is updated by multiplying the current 
value of / as 

P(A) := P(A) x / (6) 
and a histogram of the value of A is incremented as 

H(A) := H(A) + 1, (7) 

each time a configuration J with the quantity A is visited in simulations. 

These steps continue until the accumulated histogram H (A) is approximately flat. In practice, 
the histogram is regarded as flat when the desired range of A is not less than certain percent, 
say 90%, of the average value of H{A) over the range. Then, the value of the modification factor 



/ is reduced with := y/( n ) and the histogram is reset to H(A) = [191 120] . Again, the 

importance-sampling MC algorithm is performed with the new value of /. The guiding function 
P(A) is adjusted by gradually changing the modification factor / during the simulation so that 
the configurations J with the value of A are visited with equal probability. This process is 
repeated for many times until the value of / is very close to unity, e.g., / = 1.0001. 

In an early stage of simulation with a large /, the detailed balance does not hold precisely since 
the transition probability of the Markov chain and the acceptance probability for the update in 
equation ([5]) are modified by the factor / in simulations. After the factor / is sufficiently close 
to one through the iteration, the detailed balance is recovered and the stationary distribution 
of the Markov chain becomes the inverse of the guiding function, 1/P(A). Then, the desired 
distribution P(A) is estimated by 

P(A) oc P(A)H(A). (8) 



3. Griffiths singularity in randomly diluted Ising model 

The bond-diluted Ising model is defined by the Hamiltonian, 

H(S; J) = -J2 J iAS 3 -hJ2$i (9) 

(y) * 

where Si(i = 1, ■••,JV) denote Ising spins which take ±1, the interactions between the spins 
Jij are independent random variables taking J and with probability p and 1 — p, respectively 
and h represents an external field. The summation of the first term is over all the nearest 
neighbor pairs. For p close to one, a ferromagnetic phase transition occurs at a finite temperature 
T c (p). The transition temperature decreases with p decreasing and eventually it vanishes at the 
percolation threshold p c . The phase diagram of the two-dimensional bond-diluted Ising model 
is shown in figured! 
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Figure 1. Phase diagram of the two- 
dimensional bond diluted Ising ferromagnet 
in the p-T plane. The shaded region is called 
the Griffiths phase. 



According to the argument by Griffiths [14J, the free energy of the random bond Ising 
ferromagnet must be a non-analytic function of the external field h in the temperature regime 
above the transition temperature T c {p) and below T c (p = 1) of corresponding non-diluted system. 
The free energy has an essential singularity as a function of h in the zero h limit, which is refereed 
to as Griffiths singularity. Physically, the Griffiths singularity is related with the existence of 
arbitrarily large compact clusters in which most of the interactions take J. The probability of 
the existence of such large clusters is suppressed exponentially in their size and the contribution 
of such clusters to the susceptibility is proportional to the size at low temperatures. This means 
that the arbitrarily large clusters have rare but significantly important effect leading to the 
Griffiths singularity. However, the singularity of the free energy is believed to be too weak to find 
in experiments and numerical simulations. Meanwhile, such large clusters have also important 
consequences for dynamical properties because they flip from one to the other stable states very 
slowly at low temperatures. In fact, many theoretical an( ^ numerical [2HES] studies 

have found that time correlation functions in the temperature regime T c {p) < T < T c (p = 1) 
exhibit non-exponential slow relaxation which is not expected in the simple paramagnetic phase. 
Thus, the temperature regime lying between T c (p) and T c (p = 1) is called the Griffiths phase. 

In this work, we study the distribution of the (inverse) susceptibility \ at finite temperatures. 
We consider that the distribution of the susceptibility could be used as an indicator of the 




Figure 2. Schematic pictures of the distribution of the inverse susceptibility for a random 
bond Ising model, (a) At T > T c (p = 1), the distribution is bounded, (b) In the Griffiths phase 
at T c (p) < T < T c (p = I), the smallest edge of the distribution reaches the origin, (c) At the 
transition temperature T = T c (p), the statistical weight near the origin is significant and the 
average of the susceptibility diverges. 



Griffiths singularity on the analogy of the Bray-Moore argument [26] in which the density of 
states of the inverse susceptibility matrix is discussed. As is sketched in figure El the distribution 
of the inverse susceptibility P{x~ l ) is bounded at high temperatures and as the temperature 
decreases the edge of the distribution touches the origin without statistical weight. We regard 
this feature of the distribution as the Griffiths singularity. Using the MC algorithm described 
in the previous section with setting (A) = x -1 > we evaluate the distribution of the inverse 
susceptibility numerically near the origin in the Griffiths phase. 

4. Application of the importance-sampling algorithm 

We perform the importance-sampling MC simulation for sampling the interaction bonds J. In 
the inner loop of the step (ii) in section [2] we have to calculate the bulk susceptibility for a given 
J, which is expressed in disordered phase as 



For unfrustrated spin systems, cluster MC algorithms are known to be quite efficient for reducing 
critical slowing down near phase transitions \27\ 128]. Furthermore, some physical quantities can 
be calculated using an improved estimator expressed in terms of clusters, which makes statistical 
error reduced substantially. The improved estimator for the susceptibility is given by the mean 
squared cluster size in the Swendsen-Wang cluster algorithm. Using the cluster MC algorithm 
with the improved estimator [27], the susceptibility of the model is calculated with sufficient 
accuracy. 

We firstly perform the simple-sampling algorithm for 10 3 ~ 10 4 independent interaction bonds 
and obtain a first guess for the guiding function P{\ 1 )- Starting from an initial configuration 
of the interaction bonds randomly set to J with the probability p and with 1 — p, we choose a 
small percentage of the total bonds at random, and generate a new candidate for the interaction 
bonds by refreshing the chosen bonds with the same probability p. For the new candidate, we 
calculate the susceptibility by the method mentioned above and accept it with the probability 
of equation (JS]), in which a working guiding function P(x~ l ) is used and also update by the step 
(iii) in section [21 Because we are interested in the behavior of -P(x _1 ) near th origin and not in 
P{x l ) with a large value of x~ l i we put bounds to the maximum value of x" 1 allowed in the 
simulation by [x _1 ] +3cr with a being the standard deviation, which is roughly estimated in the 
preliminary simple-sampling run. In actual simulation, a new candidate with the value beyond 
the limit is rejected with probability one. 

Figure [3] shows an example of Monte Carlo trajectory of the value of x _1 for linear size 
L = 32, p = 0.6 and T/J = 1.5 and figure [H shows the histogram of the inverse susceptibility 
in this simulation. It is found that the importance-sampling MC simulation can cover a wide 
range of x ■ The re-weighted distribution by the guiding function with equation ([8]) and also 
the simple-sampling histogram for comparison are shown in figure [5] We see that the resolution 
in the simple-sampling result determined by the number of samplings is largely improved by the 
use of the importance sampling which leads to probabilities of order 10 -8 . 

We have calculated the inverse-susceptibility distributions with L = 8, 16, 32, (and 64 for 
some cases) at a few temperatures in the Griffiths and paramagnetic phases. In the paramagnetic 
phase above T c (p = 1), the distribution functions of the inverse susceptibility are scaled well 
as a function of X = (x ~ Xo)/ a with xo and cr being the average and the standard deviation, 
respectively (not shown here), and the scaling function is well described as a Gaussian. Figured] 
shows the distribution as a function of l/x -1 in the Griffiths phase with p = 0.6 and T j J = 2.0. 
In contrast to that in the paramagnetic phase, the distribution has an exponential tail for large 
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Figure 3. A Monte Carlo trajectory of the 
value of x _1 °f the two-dimensional bond- 
diluted Ising model for L = 32, p = 0.6 
and T/J = 1.5. The horizontal axis means 
renormalized MC steps which are incremented 
by one when a new value of is accepted. 
The value of x _1 as a function of the MC 
step is represented by the straight line and 
the number of rejections at the MC step is 
given by bar chart. 



Figure 4. Histogram of the inverse 
susceptibility obtained by an importance 
sampling algorithm of the two-dimensional 
bond-diluted Ising model. The parameters 
used in the simulation is the same as those 
in figure O 




Figure 5. Distribution of the inverse susceptibility of the two-dimensional bond-diluted Ising 
model with p = 0.6 for L = 32 and T/J = 1.5. The solid line represents P{\ 1 ) obtained by 
the importance-sampling MC, and the triangles are that by the simple sampling. 
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Figure 6. Distribution functions of the inverse susceptibility as a function of l/x -1 of the 
two-dimensional bond-diluted Ising model with p = 0.6 and T/J = 2.0. The system sizes are 
N = 8 2 (v)) 16 2 (<)) and 32 2 (A). The straight lines are the fitting result of the exponential 
function p(x _1 ) = B exp(— C/x^ 1 ) with B and C being fitting parameters. 



where C is a temperature dependent parameter. This behavior is also described by the 
exponential tail of a Gumbel distribution. The parameter C is obtained by fitting the function 
form of equation (jlip for large x to the data for each size and temperature. The fitting result for 
C has still size dependence. We thus extract the thermodynamic value by fitting a polynomial 
function of 1/L to the value of C. The temperature dependent of the extrapolated value of C 
for the exponential tail is shown in figure [7J In general [26J , we expect that the parameter C 
vanishes as T goes to T c (p) from above since P{x~ l ) should be a power law of x f°r T = T c (p), 
while C diverges as T goes to T c (p = 1) since P(x~ 1 ) should have a gap from the origin x" 1 = 
for T > T c (p = 1). Our result shown in figure [7J is consistent with this argument although 
further studies would be necessary to confirm the conclusion. It would be interested to see the 
precise functional form of C near T c (p) and T c (p = 1). 
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Figure 7. Temperature dependence of the 
parameter C for the left axis and \/C for the 
right. 



5. Conclusions and perspectives 

In summary, we have developed an importance-sampling MC algorithm for quenched disordered 
systems originally proposed in references [151 116] . The algorithm discussed is an efficient 
sampling method with appropriate weight for the quenched variables, which, for instance, 
correspond to the interaction bonds in random spin systems, while most of simulations have 
so far used a simple-sampling method in which the quenched variables are generated by a given 
probability distribution. Reference [TB] suggested that a learning procedure could be simplified 
by introducing a given guiding function as the weight for sampling the quenched variables. In 
this work, we explicitly formulate the learning procedure of the weight using the multicanonical 
method [I] combined with Wang-Landau recursion TP, 20] . Because the method explained in this 
paper does not assume a priori knowledge of the weight, it could be complementary to the method 
in reference [16] for studying a wide class of quenched random systems. It should be noted that 
the method can be applied to finite-temperature calculation of any distribution of observables as 
demonstrated in this work, while the applications are limited to ground-state-energy distribution 
previously. For example, the method can be in principle applicable to distribution functions of 
relaxation time and/or solving time for randomly constraint satisfaction problems. 

We have applied the algorithm to evaluate the distribution of the susceptibility in a two- 
dimensional bond-diluted Ising model which concerns the so-called Griffiths singularity. We 
have found that the distribution has an exponential tail for large value of the susceptibility and 
that the evaluated slope of the exponential tail exhibits significant temperature dependence in 
the Griffiths phase which is consistent with a simple argument [26]. It would be interesting to 
confirm the Griffiths singularity in spin-glass systems. Then, the cluster MC method used in 
the inner-loop calculation for taking an average over spin variables should be replaced by an 
extended-ensemble based MC method [3] because it does not work efficiently in frustrated spin 
systems such as spin glasses. This would require relatively extensive computational time but 
it would be helpful to use massively parallel simulations for the Wang-Landau recursion in our 
method . 

The importance-sampling MC method discussed in this paper is also regarded as a sampling 
method of rare events as originally pointed out in reference [15]. A direct application in this 
respect is to see rare samples which occur very unlikely. Figures [8] show two rare examples of 
interaction bonds found in our simulations of the bond-diluted Ising model. This makes possible 
further research on the nature of such rare samples, for example fractal properties not discussed 
here, which causes the large susceptibility and the Griffiths singularity. This method could be 
generally useful as an experimental tool for picking up the rare events and for studying their 
properties. 
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